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Abstract 

Nonequilibrium steady states in driven diffusive systems exhibit many features 
which are surprising or counterintuitive, given our experience with equilibrium sys- 
tems. We introduce the prototype model and review its unusual behavior in different 
temperature regimes, from both a simulational and analytic view point. We then 
present some recent work, focusing on the phase diagrams of driven bi-layer systems 
and two-species lattice gases. Several unresolved puzzles are posed. 



1 Introduction 



In nature, there are no true equilibrium phenomena, since all of these require 
infinite times and infinite thermal reservoirs or perfect insulations. Neverthe- 
less, for a large class of systems, it is possible to set up conditions under which 
predictions from equilibrium statistical mechanics provide excellent approxi- 
mations, as many of the inventions of the industrial revolution can attest to. 
By contrast, non-equilibrium phenomena are not only ubiquitous, but often 
elude the powers of the Boltzmann-Gibbs framework. Unfortunately, to date, 
the theoretical development of non- equilibrium statistical mechanics is at a 
stage comparable to that of its equilibrium counterpart in the days before 
Maxwell and Boltzmann. Using the intuition developed by studying equilib- 
rium statistical mechanics, we are often "surprised" by the behavior displayed 
by systems far from equilibrium, even if they appear to be in time-independent 
states. For example, the well honed arguments, based on the competition be- 
tween energy and entropy, frequently fail dramatically. In these lectures, we 
will present some explorations into the intriguing realm of non- equilibrium 
steady states, focusing only on a small class, namely, driven diffusive systems. 
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Motivated by both the theoretical interest in non-equilibrium steady states and 
the physics of fast ionic conductors [1], Katz, Lebowitz and Spohn [2] intro- 
duced a deceptively minor modification to a well-known system in equilibrium 
statistical mechanics: the Ising [3] model with nearest-neighbor interactions. 
In lattice gas language [4], the time evolution of this model can be specified 
by particles hopping to nearest vacant sites, with rates which simulate cou- 
pling to a thermal bath as well as an external field, such as gravity. If "brick 
wall" boundary conditions are imposed in the direction associated with grav- 
ity (particles reflected at the boundary, comparable to a floor or a ceiling), 
and if the rates obey detailed balance, then this system will eventually settle 
into an equilibrium state, similar to that of gas molecules in a typical room 
on earth. Although there is a local bias in the hopping rates (due to gravity), 
thermal equilibrium is established by an inhomogeneous particle density, at 
all T < oo. However, if periodic boundary conditions are imposed, then trans- 
lational invariance is completely restored so that, in the final steady state, 
the particle density is homogeneous, for all T above some finite critical T c , 
while a current will be present. Clearly, for gravity, such boundary conditions 
can be imposed only in art, as by M.C. Escher [5]. In physics, it is possible 
to establish such a situation with an electric field, e.g., by placing the d = 2 
lattice on the surface of a cylinder and applying a linearly increasing magnetic 
field down the cylinder xis. If the particles are charged, they will experience a 
local bias everywhere on the cylinder and a current will persist in the steady 
state. Echoing the physics of fast ionic conductors, we will therefore use the 
term "electric" field to describe the external drive and imagine our particles 
to be "charged", in their response to this drive. This is the prototype of a 
"driven diffusive system." Now, such a system constantly gains (loses) en- 
ergy from (to) the external field (thermal bath), so that the time independent 
state is by no means an equilibrium state, a la Boltzmann-Gibbs. Instead, it 
is a non- equilibrium steady state, with an unknown distribution in general. As 
discovered in the last decade, modifying the Ising model to include a simple 
local bias leads to a large variety of far-from-simple behavior. The scope of 
these lectures necessarily limits us to only a bird's eye view. Since a more 
extensive review has been published recently [6], we will restrict these notes 
to a brief "introduction" to this subject. Instead, we choose to include some 
developments since that review. 

For completeness, we summarize, in Section II, the lattice model introduced 
by KLS and the Langevin equation believed to capture its essence in the long- 
time large-distance regime. The next section (III) is devoted to some of the 
surprising behavior displayed by this system, at temperatures far above, near, 
and well below T c . A collection of recent developments is presented in Section 
IV. In the final section (V), we conclude with a brief summary. 
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2 The Microscopic Model and a Continuum Field Approach 

On a square lattice, with L x x L y sites and toroidal boundary conditions, a 
particle or a hole may occupy each site. A configuration is specified by the 
occupation numbers {rii}, where % is a site label and n is either 1 or 0. Oc- 
casionally, we also use spin language, defining s = 2n — 1 = ±1. To access 
the critical point, half-filled lattices are generally used in Monte Carlo simu- 
lations: J2i n i — L x L y /2. The particles are endowed with nearest- neighbor at- 
traction (ferromagnetic, in spin language), modeled in the usual way through 
the Hamiltonian: H = —4:JJ2<i,j>nirij, with J > 0. The external drive, with 
strength E and pointing in the — y direction, will bias in favor of particles hop- 
ping "downwards" . To simulate coupling to a thermal bath at temperature T, 
the Metropolis algorithm is typically used, i.e., the contents of a randomly 
chosen, nearest-neighbor, particle-hole pair are exchanged, with a probabil- 
ity min[l, e -( AH - tE )/ k BT^ Here, AH is the change in H after the exchange 
and e = (—1,0, 1), for a particle attempting to hop (against, orthogonal to, 
along) the drive. Note that these dynamic rules conserve particle number. For 
E — 0, this system will eventually settle into an equilibrium state which is 
precisely the static Ising model [3,4]. In the thermodynamic limit, it is known 
to undergo a second order phase transition at the Onsager critical tempera- 
ture T c (0) = (2.2692..) J/ks- When driven (E ^ 0), this system displays the 
same qualitative properties, i.e., there is a disordered phase for large T, fol- 
lowed by a second order transition into a phase-segregated state for low T. 
However, with more scrutiny, this superficial similarity gives way to puzzling 
surprises. In particular, T C (E) appears to be monotonically increasing with 
E, saturating at about 1.41T C (0) [7] for E > J. Why should T c (oo) > T c (0) 
be surprising? Consider the following "argument" . For very large E, hopping 
along y becomes like a random walk, in that AH is irrelevant for the rates. 
Therefore, hops along y might as well be coupled to a thermal bath at infinite 
T (apart from the bias). Indeed, recall that our system gains energy from E 
and loses it to the bath, so that any drive may be thought of as a coupling 
to a second reservoir with higher temperature. Then, it seems reasonable that 
T must be lowered to achieve ordering, since this extra reservoir pumps in a 
higher level of noise, helping to disorder the system! To date, there is neither 
a convincing argument nor a computation which predicts the correct sign of 
T C (E) — T c (0), let alone the magnitude. In the next Section, we will briefly 
review other puzzling discoveries, only a few of which are understood. 

To understand collective behavior in the long-time and large-scale limit, we 
often rely on continuum descriptions, which are hopefully universal to some 
extent. Successful examples include hydrodynamics and Landau-Ginzburg the- 
ories. Certainly, the <p 4 theory, enhanced by renormalization group techniques, 
offers excellent predictions for both the statics and the dynamics near equi- 
librium of the Ising model. Following these lines, we formulate a continuum 
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theory for the KLS model, in arbitrary dimension d. In principle, such a de- 
scription can be obtained by coarse-graining the microscopic dynamics [8] but 
we will pursue a more phenomenological approach here. Seeking a theory in 
the long-time, large wavelength limit, we first identify the slow variables of 
the theory. These are typically 

• ordering fields which experience critical slowing down near T c , and 

• any conserved densities. 

The KLS model is particularly simple since it involves a single ordering field, 
namely the local "magnetization" or excess particle density, </?(x, i), which is 
also the only conserved quantity. Here, x stands for (x±, X2, Xd-i, x d — U)- 
The last entry denotes the one-dimensional "parallel" subspace selected by 
E. Thus, we begin with a continuity equation, d t <p + Vj = 0, and postulate 
an appropriate form for the current j. In the absence of E, it simply takes 
its Model B [9] form, j(x, t) = — AV|^ + ?y(x, t), with the Landau- Ginzburg 
Hamiltonian Ti = /{^(Vy?) 2 + ^f 2 + Jfy? 4 }- As usual, r oc T — T c and u > 
0. While the first contribution to j is a deterministic term, reflecting local 
chemical potential gradients, the second term, 77, models the thermal noise. 
The noise is Gaussian distributed, with zero mean and positive second moment 
proportional to the unit matrix, i.e., < rjirjj >= 2o~5ij5(x. — x.')5(t — t'), i,j = 
l,...,d. In the presence of E, we should expect at least two modifications 
to j, namely first, an additional contribution modeling the nonvanishing 
mass transport through the system and second, the generation of anisotropies, 
since E singles out a specific direction. By virtue of the excluded volume 
constraint, the "Ohmic" current vanishes at densities 1 and 0, corresponding 
to tp — ±1. Writing £ for the coarse-grained counterpart of E, pointing along 
unit vector y, the simplest form is therefore = S(l — y? 2 ) [1 + 0(ip)]y, where 
the 0(f) corrections will turn out to be irrelevant for universal properties. 
Next, we consider possible anisotropies. Clearly, all V 2 operators should be 
split into components parallel (<9 2 ) and transverse (V 2 ) to E, accompanied 
by different coefficients. For example, the anisotropic version of the Model B 
term rVV will read (t||<9 2 + r±W\)ip, with two different "diffusion" coefficients 
for the parallel and transverse subspaces. Should we expect both of these to 
vanish as T approaches T C (E)7 Or just one of them - but which one? Recalling 
that the lowering of r, in the equilibrium system, is a consequence of the 
presence of interparticle interactions, we argue that E, stirring parallel jumps 
much more effectively than transverse ones, should counteract this effect in 
the parallel direction, having less of an impact in the transverse subspace. 
Thus, we anticipate that generically ry > r±, so that criticality, in particular, 
is marked by Tj_ vanishing at positive ry. This is borne out by the structure of 
typical ordered configurations, namely, single strips aligned with E, indicating 
that "antidiffusion" , i.e., rj_ < 0, dominates in the transverse directions below 
T c . Finally, the drive also induces anisotropies in the noise terms so that the 
second moment of their distribution should be taken as < rjirjj >= 2a.i5ij5(x — 
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x')S(t — t'). Since the transverse subspace is still fully isotropic, we define 
<j\ = ... = <Jd-i = o"||- Generically, however, we should expect er^ = <jj_ ^ cry. 
Summarizing, we write down the full Langevin equation: 



d t ip(x, t) = A{(r ± - Vi)Viv? + (r,| - d 2 )d 2 ip - 2a x <9 2 V 2 L y? + 

+^(V 2 ± + Kd 2 )ip 3 + £dip 2 }-Vr)(x,t) , (1) 

with noise correlations 

< ?7j(x, t)^(x', t') >= 2a i 5 ij 5(x - x')8(t - t') 



This equation forms the basis for the analytic study of the KLS model. Fortu- 
nately, we need not be concerned about the detailed dependence of its coeffi- 
cients on the microscopic T, J, and E. While they are in principle calculable 
within an explicit coarse-graining scheme, the properties that we already out- 
lined above suffice for our purposes. Moreover, it turns out that we can simplify 
the most general form, (1), depending on the temperature regime considered. 
We will return to this discussion in Sections 3.1 and 3.3 below. 

To conclude this Section, let us take a glance at the structure of (1). Its basic 
form is d t ip = ¥(tp,Wtp, ...) + (, with noise correlations < £(x, t)((x.', t') >= 
2<rN5(x — x!)5(t — t') and a > just a positive constant. For simplicity, we 
restrict ourselves to a scalar order parameter (p and noise "matrices" N which 
do not depend on ip (more general cases are discussed in, e.g., [10] and [11]). 
F is a functional of tp and its derivatives. Given this basic form, a steady- 
state solution P*[<f} for the configurational probability is easily found [12] 
provided ¥ is "Hamiltonian" , i.e., if it can be written as N acting on a total 
functional derivative: F = — N|^. In this case, P* is simply proportional to 
exp(— H/<j)i irrespective of the choice of N. We will refer to such a dynamics, 
in an operational sense [6], as "satisfying the fluctuation-dissipation theorem" 
(FDT) [13]. Model B falls into this category, with N = -V 2 . In contrast, 
for our driven system, F is just the expression in the {...} brackets of (1) 
and q = —Vi], so that N = —{a^d 2 + crj_Vj_). Thus, this F is clearly not 
Hamiltonian! In this fashion, our continuum theory reflects the fact that we 
are dealing with a generic non-equilibrium steady state. For completeness, 
we should point out that certain microscopically non-Hamiltonian dynamics 
may become Hamiltonian if viewed on sufficiently large length scales. The 
discussion of these subtleties [14-16], while intriguing, is beyond the scope of 
these lectures. 
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3 Surprising Singular Behavior, near and far from Criticality 



For the d — 2 Ising model in equilibrium, thermodynamic quantities are typ- 
ically analytic, except at a single point, i.e., T c . In particular, being a model 
with only short-range interactions, correlation functions are short-ranged far 
above T c , decaying exponentially with distance and controlled by a finite cor- 
relation length. Far below criticality, a half filled system in a square geometry 
will display two strips of equal width, as a result of the co-existence of a particle 
rich (dense) phase with hole-rich one. Correlations of the fluctuations within 
each phase are also short-ranged. More interestingly, the interfaces between 
the phases do represent soft degrees of freedom, being the Goldstone mode of 
a broken translational invariance [17]. One consequence is a divergent struc- 
ture factor (Fourier transform of the two-point correlation, of deviations from 
being straight). Although this behavior is singular, it is well understood and 
can be related to that of a simple random walk [18]. Only at criticality does 
the system possess non-trivial singular properties, the nature of which became 
transparent only in the 70's[19] despite Onsager's tour-de-force in 1944 [20]. 

When driven into non- equilibrium steady states, this picture changes dramat- 
ically. In particular, non-trivial singular behavior appears at all T. While the 
transition itself remains second order, its properties fall into a non-Ising class. 
Here, we will present some of these surprising features briefly, referring the 
reader to [6] for more details. 



3.1 T far above T c 



Focusing on the behavior far above criticality, let us study the two-point cor- 
relation G(r) (r =(x, y)) and its transform, the structure factor, S(k). We first 
point out that the drive clearly induces an asymmetry between x and y, i.e., 
anisotropy beyond that due to the lattice. Thus, we should not be surprised 
if the familiar Ornstein Zernike form of S (oc (1 + £ 2 k 2 )) were to become el- 
liptical rather than circular, e.g., similar to Fig. la. However, simulations [2] 
showed that there is a discontinuity singularity in S at the origin, as in Fig. lb. 
To be more precise, we may define 

R = ^ra ky ^ S(0,k y ) 
\im kx ^ S(k x ,0) 



and measure the discontinuity by R — 1. Though R does approach unity for 
T — > oo, it is about 4, for large E, even at twice the critical temperature. 
Further, it diverges as T — > T c ! We should remind the reader that, for the 
equilibrium case, S may diverge at criticality, but R is unity always. 
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(a) 



(b) 



Fig. 1. Schematic plot of structure factors, (a) Ellipses for a typical anisotropic Ising 
model in equilibrium, (b) "Butterfly" or "owl" pattern for a driven system. 

Strange as it may seem, this behavior can be understood within the con- 
text of the continuum approach, Eqn. (1). For T ^> T c , both t\\ and t_l arc 
positive. Thus, the local magnetization </?(x, t) fluctuates around a stable min- 
imum at zero, and neither the fourth-order derivatives nor the nonlinearities 
are necessary for stability. As a consequence, the behavior of the system in 
the disordered phase can be described by a much simpler, linear Langevin 
equation, namely, 

c^(x, t) = A(r ± Vl + Tj|<9 2 V - V77. (3) 



With the help of a Fourier transform, we can easily find (p(k,t) for any real- 
ization of 77, and then compute averages over the Gaussian distributed r], with 
its second moment given by (1). This forms the starting point for the discus- 
sion of the disordered phase, quantified by, e.g., equal-time structure factors 
or three-point functions [21]. 

Here, we focus on the former, S(k) =< </?(k, t)(p(— k, t) >, which can be easily 
computed: 

S(k) = ^ k * + °U k l m 
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so that R = [cr\\T±j I (cr±T\\). In order for this S to reduce to the equilibrium 
Ornstein-Zernike form, the FDT has to be invoked, which constrains R to 
unity. For the driven case, R is no longer constrained, so that a discontinuity 
develops. One consequence of such a singularity in S is that G becomes long- 
ranged, decaying as r~ 2 at large |r|. The amplitude, in addition to being oc 
(R — 1), has a dipolar angular dependence, so that an appropriate angular 
average of G is again short-ranged [6] . 

To end this discussion, let us follow G. Grinstein's argument [22] that such 
long-ranged power law decays should be expected. Starting with the dynamic 
two-point correlation G(r,t), we know that the conservation law leads to the 
autocorrelation: G(0,i) — > t~ d l 2 . In addition, being a diffusive system, we 
should expect r ~ \ft, at least far from criticality. Using naive scaling, we 
would write G(r, 0) — > |r| _d ! In other words, the generic behavior of G in a 
diffusive system is \r\~ d decay, not exponential. In this case, the equilibrium 
system is the non-generic one, in which the amplitude of this generic term is 
forced to vanish, by FDT. Our familiarity with equilibrium systems is so strong 
that, on first sight, power law decays far above T c appear quite surprising. 
Finally, note that a scaling argument of this kind cannot produce the angular 
dependence, a crucial feature of the generic singularities of driven diffusive 
systems. 



3.2 T far below T c 



Next, we turn to systems far below T c . Due to the conservation law, they 
typically display phase co-existence, so that the dominant fluctuations and 
slow modes are those associated with the interface between the phases. Again, 
for the d = 2 Ising model in equilibrium, there is a wealth of information 
on the properties of the interface [23]. In particular, being a one-dimensional 
object, the interface exhibits behavior identical to a simple random walk [18] 
at the large scales, such as widths diverging as \[L. In our case, if we study 
an interface aligned along the y axis, then, to a good approximation, we may 
specify the configuration by its position along x by the "height" function 
h(y). Known as capillary waves [24], these fluctuations have a venerable his- 
tory. Since the interface is a manifestation of broken translational invariance, 
h(y) are the soft Goldstone modes [17], so that the associated structure fac- 
tor Sh(q) = (h(q)h(—q)) diverges as q~ 2 . Of course, this property is just the 
Fourier version of the \[L divergence: Ji/ L q~ 2 dq = 0(L). Interfaces with di- 
vergent widths are also known as "rough"; only for d > 2 may interfaces in 
crystalline Ising-like models display transitions from rough to smooth phases. 

When driven, however, the interface width appears to approach a finite width 
at large distances. In particular, for E = 2 J, using lattices up to L y = 60 and 
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plotting the widths vs. LP with various values of p, we found that the curves 
did not straighten out with p as low as 0.05 [25]. As a phenomenon, roughness 
suppression is not novel, gravity being the most common example. However, 
the drive here is parallel to the interface, reminiscent of wind driving across 
a water surface, which has a destabilizing effect by contrast. Subsequently, a 
more detailed simulation study [26] of the interface with L y < 600 showed 
that S(q) ~ q~ 2 for large q, crossing over to q~ - 67 for small q (Fig. 2). Since 
Ji/Lq~ 0S7 dq = 0(1) , the small q behavior is consistent with p = 0. Though 
0.67 appears temptingly close to a simple rational: 2/3, there is no viable 
theory so far (despite two valiant attempts [27]). Finally, let us note that the 
crossover from rough to smooth occurs at the scale of q ~ E. Although no 
precise Monte Carlo analysis of this crossover has been performed, it is con- 
sistent with E having the units of 1/length. In this respect, such a length is 
comparable to the capillary length which controls the crossover in the gravi- 
tationally stabilized interface. Of course, in that case, the small q behavior is 
simply g°! 
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Fig. 2. Log-log plot of interface structure factor vs. wavevector. 
3.3 Critical Properties 

Finally, we turn to the critical region, described by t± > in (1). In contrast to 
the situation for T ^> T c , we now need fourth-order terms in Vj_ to stabilize the 
system against large-wavelength fluctuations. Similarly, the nonlinear terms 
are required to ensure a stable ordered phase below T c . However, we still 
have T|| > 0, so that fourth-order parallel gradients may safely be neglected. 
Thus, near criticality, the leading terms on the right hand side of (1) are 
(V 2 L ) 2 (f and d 2 ip, implying that parallel and transverse wave vectors, kj_ and 
k||, scale naively as |k||| ~ |kj_| 2 , i.e., parallel gradients are less relevant than 
transverse ones. More systematically, a naive dimensional analysis [28] reveals 
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that the nonlinearity associated with £ is the most relevant one, having an 
upper critical dimension d c — 5, distinct from the usual value of 4 for the Ising 
model. Dropping irrelevant terms, (1) simplifies to 

11 

c^(x, t) = \{(t ± - Vi) Viv^ + dV + ^Viv^ 3 + f V} - V ± 77 ± (x, t) (5) 



with noise term = l,...,d— 1) 

< W x > t)r) ±j (x', t') >= 2a ± 5 ij 5( y x - ar')£(t - t'). 

Here, we have rescaled t\\ to 1 and have kept the (naively irrelevant) nonlinear- 
ity associated with u, to ensure a stable theory below T c . This is the starting 
point for the analysis of universal critical behavior. 

In this regime, large fluctuations on all length scales dominate the behavior 
of the system, so that renormalization group techniques are indispensable. 
To summarize very briefly, the Langevin equation (5) is recast as a dynamic 
functional [29], followed by a renormalized perturbation expansion in e = 
d — d c [28] . The quartic coupling u must be treated as a dangerously irrelevant 
operator. Gratifyingly, the series for the critical exponents can be summed, 
so that we obtain quantitatively reliable values even in two dimensions. The 
details are quite technical [30], so that we just review the key results here. 

The discussion leading to (5) already suggests that the critical behavior of the 
driven system is distinct from its equilibrium counterpart: the upper critical 
dimension is shifted to d c = 5, and parallel and transverse wave vectors scale 
with different exponents. Anticipating renormalization, we reformulate their 
scaling as |k||| ~ |kj_| 1+A , introducing the strong anisotropy exponent A. To 
illustrate its importance, let us consider the wave vector scaling for the equilib- 
rium Ising model. For isotropic exchange interactions, coarse-graining results 
in the usual Landau-Ginzburg Hamiltonian with gradient term (Vy?) 2 . Clearly, 
this cannot lead to anything but |kj|| ~ |kj_|. The only effect of anisotropies 
in the microscopic couplings is to give rise to different amplitudes, so that 
A remains zero. We refer to this situation as weak anisotropy, in contrast to 
strong anisotropy where A ^ 0. Examples of the latter in equilibrium models 
include, e.g., Lifshitz points or structural phase transitions [31]. 

For any system with strong anisotropy, irrespective of its universality class, the 
renormalization group predicts the general scaling form of, e.g., the dynamic 
structure factor near criticality: 

5(k, t; rj.) = ^S{k\\/^, kjn, t//; T ± ///<0 (6) 
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Here, \i is just a scaling factor. Eqn. (6) can be viewed as a definition of the 
critical exponents is, z, 77 and A. The appearance of the latter is of course 
consistent with our earlier discussion. Different universality classes are dis- 
tinguished by the characteristic values of these exponents, expressed, e.g., 
through their e-expansions. For our model, all exponents, except A, take their 
mean-field values: u — z — 4, rj — while A = 1 + |. A separate analysis 
yields the order parameter exponent (3 — |. Note that all of these equalities 
are exact, i.e., all higher order terms in the e-expansion vanish! 

Care must be taken on two fronts, both associated with the presence of strong 
anisotropy, when comparing these exponents to Monte Carlo data. First, in 
the equilibrium Ising model, the same exponent rj characterizes the divergence 
of the critical structure factor near the origin in k-space, S(k, t = 0; r± = 0) ~ 
|k|-( 2 - T >) for k — > 0, and the power law decay of its Fourier transform, the two- 
point correlations, G(r, t = 0; r_|_ = 0) ~ | r |~( d ~ 2+r ?) as |r| — > 00. In contrast, 
four r^-like exponents are needed in the driven case! Fortunately, scaling laws 
relate them to rj and A. For example, we can define i]± via S(k\\ = 0,kj_,t = 
0; t_l = 0) ~ |k_i_|~( 2 ~' ?± * ) , for |kj_| — > 0. Using (6), we read off the simple 
result ?7_|_ = rj. Similarly, we introduce 77(1 through the relation G(r\\,r± = 
-(d-2+77') 

0, t = 0; t_i_ = 0) ~ r<< , for r\\ — > 00. Keeping track of A in the Fourier 

transform, we obtain the less trivial scaling law rj'^ = r? ~^^~ 3 - ) . Since most 
simulations have focused on two-dimensional systems, we set e = 3, predicting 
A = 2 and rj'^ = 2/3. Monte Carlo data for the pair correlation along the 

drive beautifully display the expected r^ 2 ^ 3 decay [32]. Numerical evidence for 
the strong anisotropy exponent is somewhat more indirect, extracted from an 
anisotropic finite size scaling analysis [7]. Convincing data collapse is obtained, 
using anisotropic systems of size L\\ x L± with fixed "aspect" ratio L||/L^ +A , 
for the theoretically predicted values of the exponents. 

It is intriguing that the signals of continuous phase transitions in or near 
equilibrium, namely, a diverging length scale, scale invariance and universal 
behavior, also mark such transitions in far-from-equilibrium scenarios. 



4 Some Recent Developments 

Beyond the topics discussed in the previous section, we may arbitrarily name 
three other "levels" of non-equilibrium steady-state systems. The first are as- 
sociated with the KLS model itself, including fascinating results on higher 
correlations [21], failure of the Cahn-Hilliard approach to coarsening dynam- 
ics [33], models with shifted periodic and/or open boundary conditions [34], 
and systems with AC or random drives [35,16]. Topics in the next "level" in- 
volve various generalizations of KLS, such as other interactions [36,37], multi- 
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layer systems, multi-species models, and systems with quenched impurities 
[38]. Further "afield" is a wide range of driven systems, e.g., surface growth, 
electrophoresis and sedimentation, granular and traffic flow, biological and ge- 
ological systems, etc. In these brief lecture notes, we present two topics at the 
"intermediate level" : a driven bi-layer system and a model with two species. 



4-1 Phase Transitions in a Bi-layer Lattice Gas 

The physical motivation for considering multi-layered systems may be traced 
to intercalated compounds [39]. The process of intercalation, where foreign 
atoms diffuse into a layered host material, is well suited for modeling by a 
driven lattice gas of several layers. In both physical systems and Monte Carlo 
simulations of a model with realistic parameters [40], finger formation has been 
observed. Both are transient phenomena, so that we might ask: are there any 
novel phenomena in the driven steady states'? On the theoretical front, there 
are two motivations. To study critical behavior of the single layer KLS model, 
it is necessary to set the overall density at 1/2, with the consequence that two 
interfaces develop below criticality. These interfaces seriously complicate the 
analysis, since their critical behavior is quite distinct from the bulk. One way 
to avoid these difficulties is to consider a bi-layer system [41], with no inter- 
layer interactions. Through particle exchange between the layers, the system 
may order into a dense and a hole-rich layer, neither of which has interfaces. 
Certainly, this is expected (and observed) in the equilibrium case, where the 
steady state is just two Ising models, decoupled except for the overall conser- 
vation law. The other is using multilayer models as "interpolations" from two- 
to three-dimensional systems [42]. Motivations aside, what the simulations 
reveal is entirely unexpected. 

In the simplest generalization, a bi-layer system is considered as a fully periodic 
L x L x 2 model. The first simulations were carried out on the case with 
zero inter-layer interactions [42]. Therefore, within each layer, all aspects are 
identical to the KLS model. Particle exchanges across the layers, unaffected 
by the drive, are updated according to the local energetics alone, so that 
these jump rates are identical to a system in equilibrium. The overall particle 
density is fixed at 1/2. In the absence of E, there would be a single second 
order transition at the Onsager temperature. Since the inter-layer coupling 
is zero, phase segregation at low temperatures will occur across the layers, 
i.e., the ordered phase is characterized by homogeneous layers of different 
densities (opposite magnetizations, in spin language). Having no interfaces, 
the free energy of the system is clearly lowest in such a state. At T = 0, one 
layer will be full and the other empty, so that this state will be labeled by FE. 

Intriguingly, when E is turned on, two transitions appear [42]! As T is lowered, 
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the homogeneous, disordered (D) phase first gives way to a state with strips 
in both layers, reminiscent of two entirely unrelated, yet aligned, single layer 
driven systems. We will refer to this state as the strip phase (S). As T is lowered 
further, a first order transition takes the S phase into the FE state. Why the 
S state should interpose between the D and FE states was not understood. 
This puzzle, together with the presence of interlayer couplings in intercalated 
compounds, motivated our study of the interacting bi-layer driven system [43]. 
Within this wider context, the presence of two transitions is no longer a total 
mystery. On the other hand, this study reveals several unexpected features, 
leading to interesting new questions. 

Our system consists of two L 2 Ising models with attractive interparticle in- 
teractions of strength unity. Arranged in a bi-layer structure, the inter-layer 
interaction is specified by J. Thus, the "internal" Hamiltonian is just H = 
—4J2 nn ' — 4:JJ2 nn ", where the first sum is over nearest neighbors within a 
given layer while n and n" differ only by the layer index. The external drive, 
E, is imposed through the jump rates, as in Section II. For simulations, we 
kept the overall density at 1/2 and chose J G [—10,10]. Note that negative 
J's are especially appropriate for intercalated materials [39,40]. The eventual 
goal of such a study is to map out the phase diagram in the T-J-E space. So 
far, we have data for only three (positive) values of E. Not expecting further 
surprises, we believe that we have uncovered the main features of the phase 
diagram. 

To begin, let us point out various features in the equilibrium case. Here, the 
J > and the J < systems can be mapped into each other by a gauge 
transformation. Thus, in the thermodynamic limit, the phase diagram in the 
T-J plane is exactly symmetric. However, for the lattice gas constrained by 
a conservation law, the low temperature states of these two systems are not 
the same, being S or FE, respectively. While the D-S or D-FE transitions are 
second order (shown as a line in Fig. 3), the S-FE line will be first order in 
nature (T axis from to 1 in Fig. 3). Due to the presence of interfaces, the 
FE domain includes the J = line. For the same reason, it "intrudes" into 
the J > half plane by 0(1/ L), for finite systems. Of course, on the T axis 
lies a pair of decoupled Ising systems, so that T c (J=0) assumes the Onsager 
value, as L — > oo. In general, T C (J) is not known exactly, though T c ( J=±oo) is 
precisely 2T C (J=0), since every configuration of the two layers can be mapped 
into one in the usual Ising model. The lesson from equilibrium is now clear: 
the S state is as generic as the FE state in this wider context. As we will see 
below, the presence of two transitions does not represent a qualitative change 
from the equilibrium phase structure. 

Turning to driven systems, we can no longer expect a symmetric phase dia- 
gram, since the drive breaks the Ising symmetry. The data for E = 25 are 
shown in Fig. 3, in which the first /second order transitions are marked as 
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Full/Empty 



Strip 



Fig. 3. Phase diagrams for the bi-layer Ising lattice gas. Solid lines are second order 
transitions in the equilibrium case. Open/solid circles are continuous/discontinous 
transitions in the driven case. 



solid/open circles. We see that the only effect of the drive is to shift the phase 
boundaries. No new phases appear while the nature of the transitions remains 
unchanged. However, there are remarkable features. One is the lowering of 
the critical temperature for large \J\. Given that T C (E) > T c (0) in the KLS 
model, it is quite unexpected that T c (\ J\ ^> 1,E ^> 1) is less than its equilib- 
rium counterpart. Perhaps more notable is the presence of a small region in 
the J < half plane, in which an S-phase is stable. Since particle-rich strips lie 
on top of each other, such a phase could not exist if either energy or entropy 
were to dominate the steady state. Concerned that its presence might be a 
finite size effect, we performed simulations using L's up to 100, with J = —0.1 
and T G [1.00, 1.20]. In all cases, the S-phase prevailed, leading us to conjec- 
ture that this region exists even in the thermodynamic limit. On the other 
hand, for lower T, the FE-phase penetrates into the J > half plane, as in 
equilibrium cases. Energetics seem to take the upper hand here. Though we 
believe that, as L — > oo, this part of the phase boundary will collapse onto 
the J = line, we have not checked the finite size effects explicitly. Of course, 
these "intrusions" into the positive/negative J region by the FE/S-phase are 
responsible for the appearance of two transitions in the studies of the J = 
model [42]. 

It would clearly be useful to develop some intuitive arguments which can "pre- 
dict" these qualitative modification of the equilibrium phase structure. Since 
the usual energy-entropy considerations appear to fail, our attempt [43] is 
based on the competition between suppression of short-range correlations [44] 
and enhancement of long-range ones [32] , as a result of driving. Let us focus on 
the two-point correlations in the disordered phase with T being lowered from 
oo and see how they are affected by the drive. On the one hand, the nearest- 
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neighbor correlations are found to be suppressed by E, consistent with the 
picture that the drive acts as an extra noise in breaking bonds. Taken alone, 
this suppression would lead to the lowering of the critical temperature, as 
pointed out in Section II. On the other hand, as we showed above, the drive 
changes significantly the large distance behavior of G, from an exponential 
to a power law. Further, the amplitude is positive (negative) for correlations 
parallel (transverse) to E. Both the positive longitudinal correlations and the 
negative transverse ones should help the process of ordering into strips parallel 
to the field. In other words, the enhanced long-range parts favor the S-phase. 
Taken alone, we expect this effect to increase T C (E). Evidently, for the single 
layer case, the latter effect "wins" . For a bi-layer system, we need to take into 
account cross-layer correlations, which are necessarily short-ranged. 

Focusing first on the J = case, where short-range effects due to cross-layer 
interactions are absent, we are led to T c (0, E)/T c (0, 0) > I. Indeed, this ratio 
is comparable to that in the single-layer case. Next, we consider systems with 
positive J. Without the drive, T c ( J, 0) is of course enhanced over T c (0,0). 
For non-vanishing drive, it is not possible to track the competition between 
the short- and long-range properties of the transverse correlations. Evidently, 
for small J, the long-range part still dominates, so that T C (J, E)/T C (J, 0) > 
1. However, for J ^> J , the presence of E effectively lowers J, since the 
latter is associated with only short-range correlations. Since a lower effective 
J naturally gives rise to a lower T c , we would "predict" that T C (J, E)/T C (J, 0) 
could decrease considerably as J increases. From Fig. 3, we see that T c ( J, E) < 
T c ( J, 0) for J > 5! The interplay of the competing effects is so subtle that 
either can dominate, in different regions of the phase diagram. Finally, we 
turn to the J < case. Here, the two effects tend to co-operate rather than 
compete, since the long-range parts favor the S-phase over the FE-phase. As 
a result, the domain of FE is smaller everywhere. In particular, note that 
the J < branch of T c ( J, E) is significantly lower than the J > branch. 
The small region of the S-phase can be similarly understood, at least at the 
qualitative level. In this picture, we may argue that the bicritical point and 
its trailing first order line should be "driven" to the J < half-plane. 

To end this subsection, we point out a few other interesting features. The 
order parameters of the S and FE phases are conserved and non-conserved, 
respectively. Based on symmetry arguments, we may expect that the critical 
properties of the two second order branches will fall into different universality 
classes. In contrast, the gauge symmetry forces static equilibrium properties 
to be identical along these two branches [45]. Since the drive breaks this sym- 
metry, we can only speculate that the D-S transitions belong to the same class 
as the single layer, KLS model [30] while the D-FE ones should remain in the 
equilibrium Ising class [14]. Since these two classes are distinct, we can also 
expect a rich crossover structure near the bi-critical point. Assuming these 
conjectures are confirmed, we can truly marvel at the novelties a driving field 
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can bring. 



4-2 Biased Diffusion of Two Species 

If the particles in the KLS model are considered "charged", it is natural to 
explore the effects of having both positive and negative charges. Defined on a 
fully periodic L x x L y square lattice, a configuration of such a "two species" 
model can be characterized by two local occupation variables, n+ and n~ , 
which equal 1 (0), if a positive or negative particle is present (absent) at site 
r=(x,y). In spin language, this corresponds to a spin-1 model, and novel 
behavior is to be expected, since we have altered the internal symmetries of 
the local order parameter. For simplicity, we will restrict ourselves to zero to- 
tal charge, Q = Y, r ( n t ~ n r) = 0- Pursuing the analogy to electric charges, 
positive (negative) charges move preferentially along (against) the drive E 
which is directed along the — y axis. As a first step, we assume that there 
are no interparticle interactions apart from an excluded volume constraint. 
Thus, our model can be viewed as the high-temperature, large E limit of 
a more complicated interacting system. To model biased diffusion, a parti- 
cle with charge q — ±1 jumps onto a nearest-neighbor empty site accord- 
ing to min[l, e~ qESy ], where 5y = 0, ±1 is the change in the ^-coordinate 
of the particle. To allow charge exchange between neighboring particles, two 
nearest-neighbor sites carrying opposite charge may exchange their content 
with probability 7min[l, e~ qESy ]. Now, 5y is the change in the ^-coordinate of 
the positive particle. The parameter 7 sets the ratio of the characteristic time 
scales controlling, respectively, charge exchange and diffusion. 

Physical motivations for this model come from various directions: fast ionic 
conductors with several species of mobile charges [1], electric breakdown of 
water-in-oil microemulsions [46], and gel electrophoresis of charged polymers 
[47]. It can also be interpreted as a simple model for some traffic or granular 
flows [48,49]. 

Summarizing our simulation data, we first map out the phase diagram of the 
model, in the space spanned by E, the total mass density m = J2r{ n t + 
n~) / (L x L y ) , and 7 [50]. Focusing on small 7, the diffusive dynamics is limited 
by the excluded volume constraint. For sufficiently small E and rh, typical 
configurations are disordered, characterized by homogeneous charge and mass 
densities and fairly large currents. Nevertheless, this phase is highly nontrivial, 
supporting anomalous two-point correlations reminiscent of the KLS model: In 
generic directions, the familiar r~ d decay prevails, with the remarkable excep- 
tion of the field direction, where a novel r^ d+1 ^ 2 dominates [51]. With increas- 
ing E or fh, the tendency of the particles to impede one another becomes more 
pronounced, until a phase transition into a spatially inhomogeneous "ordered" 
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phase occurs. Beyond the transition line E c ( m,7), typical configurations are 
charge-segregated, consisting of a strip of mostly positive charges "floating" 
on a similar strip of negative charges, surrounded by a background of holes. 

To explore these transitions more quantitatively, we need to define a suitable 
order parameter. Taking the Fourier transforms of the local hole and charge 
densities, (j)(k±,k\\) = j^j- X) r [l ~~ ( n t + n r)] exp{ik±x + ik\\y} and ip(k±, k\\) = 

7777- ErK ~~ n r] Gxp{ikj_x + ik\\y}, we select the components $ = |0(O, |f-)| 
and ^ = |^(0, |^)| since these are most sensitive to ordering into a single 
transverse strip. Naturally, we choose their averages < $ > and < ^ > as our 
order parameters. To distinguish first from second order transitions, we also 
measure their fluctuations and histograms. The resulting phase diagram is 
shown in Fig. 4. Typically, hysteresis loops in < $ > , < ^ > and the average 
current are observed only for fh < m ( r y), indicating that the transitions are 
first order in this regime. On the other hand, fluctuations develop sharp peaks 
for larger mass densities, signalling continuous transitions. Larger values of 7 
favor the disordered phase, so that the whole transition line shifts to higher 
E, the region between the spinodals narrows, and m {j) moves to higher mass 
densities. Thus, we observe a surface of first order transitions, separated from 
a surface of continuous ones by a line of multicritical points: E c (m ( , y), 7). 




Fig. 4. Phase diagram for the driven two-species model. The system size is 30 x 30. 
The filled circles mark the lines of continuous transitions, while the open circles 
denote the spinodal lines associated with first order transitions. The inset shows a 
typical ordered configuration at fh = 0.40, E = 3.00 and 7 = 0.02 with the open 
(filled) circles representing positive (negative) charges. 
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While it is not surprising that the disordered phase is stable for large E pro- 
vided the mass density is sufficiently small, it may appear rather counterintu- 
itive that it should also dominate the m < 1 region. However, setting m = 1 
eliminates every hole in the system so that the dynamics is carried entirely 
by the charge exchange mechanism. Relabelling positive charges as "particles" 
and negative ones as "holes" , the model becomes equivalent to a noninteract- 
ing (i.e., J = 0) KLS model. Also termed the asymmetric simple exclusion 
process (ASEP) in the literature, its steady-state solution is exactly known 
[52] to be homogeneous. Thus, for all 7 > 0, our two-species model remains 
disordered along the entire rh = 1 line. For sufficiently large 7, the data in- 
dicate quite clearly that a finite region of disorder persists, for any E, just 
below complete filling. For smaller 7, however, it is less obvious whether such 
a region exists in the thermodynamic limit, since even a single hole can suf- 
fice to induce spatial inhomogeneities in a finite system: Acting as a catalyst 
for the charge segregation process, the hole creates a strip of predominantly 
positive charge, separated by a rather sharp interface from a similar, nega- 
tively charged domain, located "downstream". The charge exchange mecha- 
nism tends to remix the charges but has little effect for small 7. Eventually, 
the hole ends up trapped in the interfacial region! Finally, we turn to larger 
values of 7, specifically, 7 > 0.62: here, the charge exchange mechanism dom- 
inates over the excluded volume constraint and suppresses the ordered phase 
completely. 

So far, our discussion has mostly drawn upon Monte Carlo results. It is grat- 
ifying, however, that the same qualitative picture emerges from a mean-field 
theory, based upon a set of equations of motion for the local densities. We 
briefly summarize the analytic route. Since the numbers of both positive 
and negative charges are conserved, we begin with a continuity equation for 
the coarse-grained densities p ± (r,t): d t p ± + Vj ± = 0. For simplicity, we fo- 
cus on the case 7 = 0, i.e., particle-hole exchanges, first. As for the KLS 
model, ^ can be written as the sum of a diffusive piece and an Ohmic term, 
j ± = A^ — Vj^p ^ ± Sy}. The density-dependent mobility X^ 1 must vanish 

with both /r* 1 and the local hole density, ip = 1 — (p + +p~), i.e., X^ 1 = p ± <p. The 
"Hamiltonian" 7i = {p + In p + + p~ In p~ + (p In (/?} is just the mixing entropy 
associated with distributing L x L y p + positive and L x L y p~ negative charges 
over a lattice of L x L y sites. Note that the functional derivative is taken 
at fixed p T since we are focusing on particle-hole exchanges here. To model 
the charge exchange mechanism, we simply add a similar term to j ± , namely, 
j± = A{-V^r ± Sy} + 7 A'{-V^ ± Sy} where A' = p + p~ vanishes with 

both p + and p~ , and the functional derivative is taken at fixed hole density. 
Expressing the resulting equations in the more convenient variables ip and 
charge density ip = p + — p~ , we obtain 



d t (f = V{Vy? + Stpipy} 
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d t tp = V{7V^ + (1 - 7) [y?W> - ^Vy?] -£(p(l- cp)y 

-Im-vf-^m 



(7) 



These equations have to be supplemented by periodic boundary conditions 
and constraints on total mass and charge, (1 — m)L x L y = J d d r(p(r,t) and 
= / d d rip(r,t). Also, to be precise, all Laplacians V 2 should be given the 
appropriate anisotropic interpretation. 

Using (7) as our starting point, we seek to recapture the major features of the 
phase diagram. First, the presence of two phases, uniform vs. spatially struc- 
tured, is reflected in the existence of both homogeneous and inhomogeneous 
steady-state solutions of (7). The existence of the former is evident, due to 
the conservation law. Anticipating homogeneity in the transverse directions, 
we seek the latter in the form <p(y), ip{y)- Integrating (7) once, the stationary 
mass and charge currents appear as natural integration constants. Since the 
numbers of positive and negative charges are equal, the former must be zero, 
leaving us with the latter: J. The first equation allows us to eliminate ip in 
favor of ip, and, with u = yT + as the new variable, we can reduce the 

second equation to potential form, u" = — dV ^ . The potential V(u), given 
explicitly in [50], exhibits a minimum for a range of J's, so that inhomoge- 
neous, periodic solutions of (7) exist. We note in passing that these may even 
be found explicitly, provided 7 = [53]. Otherwise, numerical integration is 
of course always possible, yielding rather impressive agreement with simula- 
tion profiles. Interestingly, our equations of motion predict that £ enters only 
through the scaling variable £L y . If plotted accordingly, our simulation data 
confirm this scaling by collapsing rather convincingly, at least within their 
error bars. 

However, the mere existence of two types of steady-state solutions is not suf- 
ficient to provide evidence for a phase transition. Therefore, we seek insta- 
bilities of these solutions, as the control parameters £, fh and 7 are var- 
ied. Performing a linear stability analysis, we find that a homogeneous so- 
lution with mas s fh becomes un stable if the drive £ exceeds a threshold value 
£#(771,7) = j~ \J (i_m)[(2-7)m-i] "- The most relevant perturbation is associated 
with the smallest wave vector in the parallel direction, (0,2ir/L y ). A simple 
analysis of £#(771,7) shows that this instability can only occur within the in- 
terval 2~ < to < 1, so that, in particular, the homogeneous phase is always 
stable for 7 ^ 1. Clearly, we cannot identify this mean-field stability boundary 
with the true transition line: we have not considered the locus of instabilities 
associated with inhomogeneous solutions, and fluctuations have been neglected 
throughout. However, it is quite remarkable how well it mirrors the qualitative 
shape of the phase diagram. 

Finally, let us consider the order of the transitions. In principle, two routes 
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can be pursued here. One of these, namely, the computation of the stability 
boundary of the inhomogeneous phase, Ej(rh, , ~f), is rather subtle, involving 
three Goldstone modes [54]. Once £/(m, 7) is known, values (m, 7) for which 
£1 and Eh coincide can be identified as loci of continuous transitions; other- 
wise, Eh and £/ mark the "spinodals" , near a first order transition, where the 
homogeneous/inhomogeneous phases become linearly unstable. Alternatively, 
the adiabatic elimination of the fast modes results in an effective equation 
of motion for the slow ones which can then be analyzed. Since the technical 
details of this approach can be found in [50], we need only focus on the result 
which combines the expected with the surprising. Letting M(q±,t) denote the 
complex amplitude of the unique slow mode, associated with the band of wave 
vectors (q±,27i/L y ), its equation of motion takes the Ginzburg-Landau form: 
d t M = -{(t + q\)M + gM\M\ 2 + 0(M 5 )}. Here, r is the soft eigenvalue, 
which vanishes on the stability boundary and g = g(E,m,^,L y ) is a rather 
complicated function. As in standard Landau theory, the sign of g determines 
the order of the transition: if positive, the transition is continuous while neg- 
ative g signals a first order one. Since r ~ 0, we can set E = Eh and discuss g 
on the stability boundary itself: being positive for m < 1, it has a unique zero 
at a critical m ( , y), below which it becomes negative. Since m G (7) increases 
with 7, we recover the qualitative behavior of the multicritical line observed 
in the simulations. The surprising aspect of this equation of motion, however, 
resides in the fact that M is 0(2)-symmetric and q± spans just a single spatial 
dimension. Given these symmetries, the Mermin- Wagner theorem [55] should 
forbid the existence of long-range order! One might hope that a careful analy- 
sis of finite-size effects in this system would contribute to the disentanglement 
of this puzzle. 

We conclude with two comments. First, we noted earlier that the case rh = 1 
is exactly soluble. There are, in fact, two other surfaces in the phase dia- 
gram, corresponding to 7 = 1 and 7 = 2, for which certain distributions 
are exactly known. Setting 7 = 1 implies that the rates for particle-hole and 
particle-particle exchanges become equal, so that any given particle, e.g., a 
positive charge, cannot distinguish between negative charges and holes. Thus, 
it experiences biased diffusion, equivalent to the non-interacting KLS model. 
Accordingly, the marginal steady-state distribution of the occupation numbers 
of one species is uniform, i.e., P [{nf}} = J2{ n f} P K n r ' ■> n r}] ^ 1> an d observ- 
ables pertaining to a single species are trivial. For example, the two-point 
correlation functions for equal charges, < n^n^ > and < n^n^ >, vanish for 
r 7^ 0. The full distribution, P [{n^,n~ }], however, is nontrivial, so that, e.g., 
the two-point function for opposite charges, < n+ng >, remains long-ranged 
[51]. A completely random system, marked by a uniform P [{n+,n~}] results 
only if 7 = 2. Finally, we note that the one-dimensional version of our model 
is exactly soluble by matrix methods [56]. 

Second, it is natural to wonder about the consequences of having non-zero 
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charge. This problem has only been investigated for 7 = [57], but the find- 
ings are quite remarkable. The system still orders into a charge-segregated 
strip, while supporting a nonvanishing mass current, reflected in an overall 
drift. If, e.g., positive charges outnumber negative ones, one might expect that 
the whole strip would drift in the field direction - in analogy with American 
football, where the team with fewer players tends to lose ground. In contrast, 
the strip wanders against the field, following the preferred direction of the 
minority charge! 

We should add that this model possesses several other intriguing properties, 
e.g., stable configurations with nontrivial winding number ("barber poles") 
[58] or multiple- valued currents in the mean- field description [53]. To sum- 
marize, the remarkable richness of this deceptively simple system is clearly 
amazing. 



5 Summary and Outlook 

We have presented, within the limited scope of these lecture notes, a brief 
introduction to the statistical mechanics of driven diffusive systems and some 
recent developments in this ever expanding field. Focusing only on the proto- 
type model and two of the simplest generalizations, we pointed out a multitude 
of surprises, when we base our expectations on the experience with equilibrium 
systems. While some of these, e.g., the generic presence of singular correla- 
tions, are rather well understood, others, such as the nature of ordering in 
two-species models, remain unresolved. Of course, the holy grail of this whole 
field, namely, the fundamental understanding and theoretical classification of 
non-equilibrium steady states, still beckons at the distant horizon. 
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